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Abstract 

Recent analyses have shown that the inclusion of electroweak corrections can alter significantly 
the energy spectra of Standard Model particles originated from dark matter annihilations. 
We investigate the important situation where the radiation of electroweak gauge bosons has 
a substantial influence: a Majorana dark matter particle annihilating into two light fermions. 
This process is in p-wave and hence suppressed by the small value of the relative velocity 
of the annihilating particles. The inclusion of electroweak radiation eludes this suppression 
and opens up a potentially sizeable s-wave contribution to the annihilation cross section. We 
study this effect in detail and explore its impact on the fluxes of stable particles resulting 
from the dark matter annihilations, which are relevant for dark matter indirect searches. We 
also discuss the effective field theory approach, pointing out that the opening of the s-wave 
is missed at the level of dimension-six operators and only encoded by higher orders. 



1 Introduction 



In a recent paper [1] some of us have pointed out that the energy spectra of the Standard Model 
(SM) particles originating from Dark Matter (DM) annihilation/decays can be significantly 
affected by Electro Weak (EW) corrections, if the mass M of the DM particles is larger than 
the EW scale (here set to be the mass of the W boson m w ). The emission of EW gauge 
bosons from the final state of the annihilation/decay process is enhanced by single logarithms 
In M 2 /m 2 v in the collinear region and by double logarithms In 2 M 2 /m 2 v when both collinear 
and infrared singularities are present, and implies that all stable particles of the SM appear in 
the final spectrum, independently of the primary annihilation/decay channel. The inclusion 
of EW corrections seems therefore an essential ingredient in indirect searches for DM. The 
impact of EW corrections is particularly relevant in two situations: (1) when one is interested 
in the energy region of the final fluxes (after propagation from the source) which corresponds 
to the low-energy tail of the spectrum, populated by the decay products of the additional 
gauge bosons; (2) when some of the stable species appear only if EW corrections are taken 
into account, for instance antiprotons (from W or Z decays) in an otherwise purely leptonic 
annihilation. 

One basic assumption made in Ref. [1] was that the tree-level 2 — > 2 annihilation cross 
section of the DM particles into SM states was dominant over the 2—7-3 cross section with 
soft gauge boson emission from the external legs, and the latter was factorized with respect 
to the former. 

While this assumption is certainly reasonable and commonly made (for instance if the DM 
is a heavy Dirac fermion singlet under the SM gauge group), there are well-motivated cases in 
which it is questionable. Consider for instance a DM particle x which is a Majorana fermion 
and a SM singlet. The cross section of XX ~~ >" J I *> where DM annihilates into SM fermions 
of mass rrif, consists of a velocity-independent (s-wave) and a velocity-dependent (p-wave) 
contribution 

va = a + bv 2 + 0(v 4 ) , (1.1) 

where v ~ 10 -3 is the relative velocity (in units of c) of the DM states in our Galaxy. By 
helicity arguments, a oc (rrif/M) 2 and hence very suppressed for light final state fermions 
(e.g. leptons), while the p-wave term is suppressed by v 2 . In this case, it is not guaranteed 
that the 2-body annihilation cross section is quantitatively larger than the one with EW 
corrections. Indeed, the latter ones may open a sizeable s-wave contribution and elude the 
suppressions. 

The scope of this paper is therefore to generalize the results of Ref. [1] to the interesting 
case in which the 2-body annihilation cross-section is not automatically larger than the one 
with soft gauge boson emission. The same kind of effect has been considered in the past with 
respect to photon radiation in Refs. [2-4]; this work is partly, but not only, an extension of 
those analyses to include also W, Z gauge bosons. On the other hand, an approach similar to 
ours is the one carried out in Ref. [5], but our final results disagree with the ones published 
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there 1 . Other works which considered at various levels the impact of EW corrections on 
DM annihilation or cosmic ray physics include [6]. For somewhat related work on 3-body 
annihilations below threshold see also Ref. [7]. 

We shall show in detail that, whenever the dark matter annihilation occurs by exchange 
of a heavy intermediate state, at the lowest order in the expansion in inverse powers of this 
heavy mass, final state radiation is not sufficient to remove the helicity suppression, while it 
is efficiently removed at higher orders by processes involving the emission both from external 
legs and from virtual internal propagators. Although subleading in terms of powers of the 
heavy mass, these contributions do not pay the velocity suppression and actually can be 
dominant. 

This allows us to raise an important and cautionary remark concerning the use of the 
effective field theory approach to describe DM interactions [8]. If the interactions of DM 
with SM fermions are described by effective four-fermion dimension-six operators, then the 
emission of soft gauge bosons can only take place at the lowest order from the external 
legs and the corresponding cross section remains helicity suppressed. In other words, the 
effect of opening up a large s-wave annihilation channel is missed at the level of dimension- 
six operators and then one may (incorrectly) conclude that the whole cross section is still 
suppressed. Instead, as we shall point out, the diagrams leading to s-wave contributions 
correspond to operators with dimension higher than six, whose quantitative relevance for the 
cross section can be comparable or larger than that due to dimension-six operators, despite 
the larger dimensionality. 

The plan of the paper is as follows. In Section 2 we describe the simple model we shall use 
throughout the paper. Then, Section 3 introduces preliminary considerations about 2-body 
annihilations, especially about helicity-suppression, setting the ground for the subsequent 
discussion. Section 4 contains the calculations and the results for the annihilation cross 
section with the inclusion of EW bremsstrahlung (including the Ward Identities check and 
the remarks on the effective field theory approach). Our results for the 3-body cross sections 
are in disagreement with those of Ref. [5], while we find a perfect agreement with those in 
Ref. [2] concerning the radiation of one photon. With these analytical results at hand, we are 
then ready to study in detail the impact of the opening of the s-wave on the fluxes of stable 
SM particles. To this end, we carry out a numerical analysis. In Section 5, we derive the 
energy spectra at production while in Section 6 these spectra are then propagated to give the 
fluxes of particles at detection. Concluding remarks are collected in Section 7. 

2 The model 

In this section we present the (toy) model we shall use in the paper to describe the relevance 
of the EW corrections in DM annihilations. Let us add to the particle content of the SM a 
Majorana spinor \ with mass M x , singlet under the SM gauge group and playing the role of 

1 The disagreement originates from the use of incorrect Fierz identities, thus invalidating the calculation of 
the cross sections, as we have been informed by the authors of Ref. [5] in a private communication. 
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DM, and a scalar SU (2)-doublet S, with mass M s > M x 

x = x c «=(5)- (21) 

The field S provides the interactions of the DM with the generic fermion of the SM, described 
by the left-handed doublet L = (/i,/^)- In fact, the total Lagrangian of the model is (see 
also Ref. [9]) 

£ = £sM + £ x + £-S + Ant , (2.2) 

where to the Standard Model Lagrangian £sm we added 

£ x = \x(i9-M x ) X , (2.3) 
C s = (D,S)\D^S)-M 2 S S^S, (2.4) 
Ant = y L x(Ltcr 2 S) + h.c. = y L (xP L f2V + -xP L fiV°)+^., (2.5) 

where the 4-component notation has been used and where contractions on SU(2) indices is 
defined as (Lio^S 1 ) = Li{ia2)ijSj. Moreover, we shall adopt the convention for projectors: 
Pr,l — (1 i 7 5 )/2. The stability of the DM can be achieved e.g. by endowing \ an d S with 
odd parity under a Z 2 symmetry, while the rest of the SM spectrum is even. 

The model is manifestly gauge invariant, and is the same of Ref. [5], which will allow us a 
direct comparison between their results and ours. A reader expert in supersymmetry would 
recognize the same interactions of a Bino with fermions and their supersymmetric scalar 
partners. 

We shall restrict our attention to the massless limit — mf 2 =0. While reasonable for 
leptons and light quarks, this approximation may not be good for heavy quarks. For instance, 
the DM annihilation into tt, if kinematically allowed, would proceed through s-wave with a 
contribution proportional to (m t /M x ) 2 , which can be large already without EW corrections. 
However, the generalization of our calculations to non-zero fermion masses is beyond the 
scope of this paper. 

As anticipated in the Introduction, one of the main goals of the present paper is to 
show that the inclusion of higher-order processes with emission of soft weak gauge bosons 
evades the helicity suppression and turns on an unsuppressed s-wave contribution to the DM 
annihilation cross section. Before turning to the details of the calculation of the 2 — > 3 
scatterings in Section 4, let us first review some standard material about 2 — > 2 annihilations, 
with particular emphasis on the role of helicity suppression. This will serve to set the notation 
and to highlight the main points for later use. The analogous results for the amplitudes and 
the cross sections of the 2-body and 3-body processes in the case of Dirac DM are reported 
in Appendix A. 
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Figure 1: Feynman diagrams for the tree-level annihilation in Eq. (3.1) together with its effective 
contraction in the limit M$ S> M x . 



3 Two-body annihilation into fermions and the helicity suppression 

Let us consider the annihilation of the DM Majorana fermion into a pair of massless left- 
handed fermions (see Figure 1) 

x(h)x(h)^fLMfLM- (3.1) 

The cross section admits the usual expansion in powers of the relative velocity v of the initial 
DM particles 

vcr = a + bv 2 + 0(v 4 ) , (3.2) 
where the coefficients a, b corresponding to s- and p-waves, respectively, are given by 

a = ' 6= 487(TT7Fm|' (3 ' 3) 

where we have defined 

A 

This result shows the well-known fact that the first non-zero contribution to the tree-level 
cross section for the Majorana DM is velocity dependent, and hence suppressed. 

Let us try to understand this fact in simple terms. At the level of Feynman diagrams 
the Majorana nature of the DM implies the presence of two crossed channels, t and u (for 
Dirac DM, only the t-channel contribution would be present). The amplitudes are of the 
form (x-Pl/)(/-Prx), which becomes (xiaPRX)(f / y a PLf) after chiral Fierz transformation. 
The total tree- level amplitude for the process in Eq. (3.1) is given by 
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■Mo = t -^-[Uf(p 1 ) l a P L Vf{p2)} 



«x( fc l) + o W2)7 75«x( fc l) 



(3.5) 

where we have defined the quantities 
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which satisfy the property D a — D i2 — 2p i -(k 1 — k 2 )D n D i2 . Notice that the momenta of the 
incoming DM particles are such that k^ — k!^ ~ 0(v)M x . We thus obtain that 

D lX -D a ~0{v)M 2 x D lX D l2 . (3.7) 

In the matrix element in Eq. (3.5), the first term represents a vector current while the second 
term is an axial- vector current. Let us analyze the velocity factors present in each of them, 
in the non-relativistic limit v < 1. The vector current is multiplied by a factor proportional 
to v due to Eq. (3.7). For the axial current of Eq. (3.5), using the Gordon identities we have 

v x {k 2 )l a l 5 u x {k 1 ) = — l ^^v x {k 2 ) l5 u x {k l ) - -^-v x (k 2 )o a \k xfi - k 2 p) l5 u x {h) (3.8) 

X X 

The vector (ki + k 2 ) a = {pi+p 2 ) a in the first term saturates the current Uf'jaPLVf in Eq. (3.5) 
and gives rise to terms proportional to the fermion mass, which are zero in our computation. 
The second term gives again an 0{v ) contribution. We thus recovered the well-known fact 
that for Majorana fermions the scattering amplitude is proportional to the first power of the 
relative velocity of the incoming particles. Notice that for Dirac DM Eq. (3.5) would not 
contain the D 2 i terms, as only the t-channel contributes to the amplitude, and the vector 
current thus gives rise to an unsuppressed s-wave term in the cross section (see App. A for 
details). 

Another interesting limit to analyze is the large scalar mass regime r ^> 1 for which 
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1 + 



™< l D <-J ~ D '2 - °[^)^p- (3-9) 

X 



In this case, the amplitude for DM Majorana annihilation into (massless) fermions at leading 
order in v and 1/r is given by Eq. (3.5), where the first term in square brackets is 0(v/r 2 )/M x , 
which is subdominant with respect to the second one, of order [0(v/r)/M x ][v x a a3r y5U x ]; thus, 
the tree-level cross section will approximately be given by 

Mxx^ff)~~. (3.10) 
4 Three-body DM Annihilation 

Let us now turn to analyze the case of interest, namely the emission of EW gauge bosons 
in DM annihilations. First, we are going to manipulate the matrix element and discuss its 
velocity dependence. Then, we deal with the kinematical constraints of the 3-body phase 
space and arrive at the results for the cross section. Finally, we re-interpret our findings in 
the language of effective field theory and make some remarks about its use. 

4.1 Matrix element and velocity dependence 

Let us discuss for definiteness the 3-body process with the emission of a Z boson 

x(ki)x(k 2 ) ->■ f L (p 2 )f L (p l )Z(k) . (4.1) 
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Figure 2: Feynman diagrams for the 3-body process in Eq. (4-1)- 



The corresponding Feynman diagrams are depicted in Figure 2. Of course, we shall include 
also the emission of W gauge bosons in the final results. The amplitude can be written as 

iM . e * = tg\y L \ 2 {l-2s 2 w ) _ + ^ _ M ^ + ^ _ M ^ ^ (4 2) 

where we have denoted sw = sm9w and cw = cos 9w, being the Weinberg angle. Follow- 
ing e.g. Ref. [3], we shall call "FSR" (final state radiation) the processes where a gauge boson 
is radiated from an external leg, while we refer to the emission from internal virtual particles 
as "VIB" (virtual internal bremsstrahlung) . Thus, the A and C terms are of FSR type, while 
the B terms are VIB. 

In order to show how the Fierz transformation works, let us analyze in more detail the 
amplitude A4a, for massless outgoing fermions 

_ 2[u f { Pl )f{k)P L {f i + ffKfo)] [v x (k 2 )P L v f (p 2 )] 
A ~ (2 Pl ■ k + m|)(M2(l - r) - 2p 2 ■ k 2 ) ' 1 ' ] 

where each fermionic current is composed by a Dirac and a Majorana spinor, for / and x 
respectively. Applying the Fierz transformation 

(Pr) V {Pl)u = \{P R rUPLl») kv (4.4) 
we can perform the rearrangement 

v x (k 2 ) i (P L ) ij v f (p 2 ) j [u f (p 1 )Y(p 1 + %)]k{Pii)kiu x {kx)i 
= ^K(A; 2 )P L7 %(A; 1 )][« / (p 1 )7 p ^ 1 + flPjiwCpa)] • (4-5) 
so that Eq. (4.3) becomes 

[u f (Pl)r(k)(f 1 + ^)PB^Vf(P2)} [v x (k 2 )P L rU x (h)] 

A (2 Pl ■ k + m|)(M x 2 (l - r) - 2p 2 ■ k 2 ) ' 1 ' ] 
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With the same technique, the terms in the total amplitude (4.2) read 

KA ..exc u f f{k){f +$)P R ^v f (D 22 -D 2l _ D 22 + D 21 _ 
M A -M A = 2pi . fc + m| {— v xl »u x + — W5% 

(4.7) 

M B -MT = (-l)[u f P R Yv f ][(k 1 -k 2 - Pl +p 2 )-e*(k)v x P Llll u x D 11 D 22 

- (k 2 -h- Pl + pa) • e*(k) v xlil P L u x D 21 D 12 } , (4.8) 

Mc-M c = 2p2 . k + ml { 2 + — v ^ u * 

(4.9) 

Let us now discuss the limit ii < 1 and r ^> 1, in analogy with the previous section for the 
2-body process. The coefficient of the DM vector current in the full amplitudes of the kind 
A, B and C is always 0{v) as it happens for the 2-body case, due to Eq. (3.7); in particular, 
in the large r limit it is proportional to 0(v/r 2 )/M 2 (as in Eq. (3.9)). Instead, for the axial- 
vector current, the crucial point is that the mass cancellation of order 0(rrif) does not occur 
anymore. Indeed, from the Gordon identities we recover Eq. (3.8), where the second term 
turns out to be proportional to v, while in the first term the 4- vector saturating the fermionic 
currents is now (pi + p 2 + fc) M , which does not trigger anymore the chiral identity, and leaves 
non-zero terms even for vanishing to/. Thus, the terms of the amplitudes containing v x ^u x 
read 

(Ma _ mt + m c - mt% m = M w^* 1(g22+g2i) : ( ° 11+gl2) 

X v v ' 

\ r 2 I M 2 

(4.10) 



, , * m (JJ n JJ 22 + D l2 D 21 ) 
(p 2 - Pi) • e (fc) + 



(fci -k 2 )-e (k) 



o(4) A 

V r 2 I 

(4.11) 



\r A / Mi 



where we have highlighted the behavior of each term with v and 1/r. Notice that now there 
appear terms without v dependence. Indeed, in the limit v — and to leading order in 
1/r <C 1, the full amplitude is given by 

^Uo = %^^[(%^/)(Pi-P2)-(^i + M-(%fe/) (P2-Pi)-e*], (4.12) 

where the first term comes from FSR while the second originates from VIB, and they are both 
of order 0{l/r 2 ). Schematically, the various contributions to the amplitude can be organized 



7 



as follows 
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(4.13) 



At this point we can learn an important lesson (see also Ref. [2]): the opening of the s-wave 
originates from diagrams of both FSR and VIB type, at 0(l/r 2 ) in the amplitude; limiting 
the expansion up to 0(1/ r) in the amplitude would cause the process to stay in the p-wave. 

An order-of-magnitude estimate for the 3-body cross section, showing the leading depen- 
dence on the expansion parameters, can be obtained straightforwardly 



v(r(xx -»■ ffZ) 



aw 
Ml 



O 



+ o 



o 



(4.14) 



where the weak coupling aw = Q 2 / (47r) for the gauge boson emission has been restored. The 
estimates in Eqs. (3.10) and (4.14) allow to gather an understanding in simple terms of the 
situation we are studying. While the 2-body annihilation cross section behaves like v 2 /r 2 , 
the 3-body FSR and VIB diagrams give rise to both s-wave and p-wave terms. The p-wave 
from 3-body processes cannot compete with the 2-body cross section because of the extra 
aw factor; however the s-wave from 2 — ?- 3 annihilation, free from the v 2 suppression, can 
overcome the 2—^2 cross section if r is not too large. In the next subsection we shall give a 
more precise estimate based on the analytical results. 

Because of the importance of this point, and being the distinction between FSR and VIB 
not able to disentangle clearly the s-wave contribution from the p-wave one, let us introduce 
now a specific notation. Having in mind an expansion in powers of 1/r in the amplitude - as 
sketched in Eq. (4.13) - we shall call "leading order" (LO) the lowest order term 0(l/r) in this 
expansion, which originates from lowest order FSR-type diagrams. This is the order at which 
Refs. [1,5] work. As shown above, in the LO approximation the annihilation cross section 
proceeds through p-wave. Only higher order terms are able to remove the helicity suppression. 
We shall further elaborate on this expansion as an operator expansion in Sect. 4.3. 

As a check of the results of this subsection, one can use the Ward Identities for EW 
SM gauge bosons k^M.^ ~ 0, for rrif ~ 0, where is the amplitude computed for the 
longitudinal mode of the Z. By direct calculation one obtains 



k,(M» A -MT xc ) 

k,(M» c -Mr c ) 
k,(M» B - M^ c ) 



(ufy a P L v f ) 
-(ufj a P L v f 



Doo-D 



22 — 



v xl u x ) + 



Doo + D 



22 T -1^21 /_ 



"K7 Ux) + o K7 15U X ) 



2 v *■' 2 
= -(u fla P L v f ) y\D 22 - D n - (D 21 - D 12 )](v xl a u x ) + 

+ [^22 - D n + (D 3 i - D 12 )] (v x7 a -f 5 u x ) 



(4.15) 
4.16) 

(4.17) 



whose vanishing sum confirms the Ward Identity 



2 In the B diagrams, we used the trick (ki—k2+P2—Pi) , k = D^—D^ and {k2—ki+p 2 —pi)-k — D l2 — -D 2 i • 



S 



It is interesting to see the level of cancellation in the large Ms limit using the properties 
of Eq. (3.6). We see that, up to order 0(1/ r), M. A and M. c cancel each other so that, at this 
order, we can say that their sum is gauge invariant; if we want to keep corrections of order 
0(l/r 2 ) or higher, the full sum of A, C and B diagrams have to be considered in order to 
have a consistent evaluation 

k,(M» A - MT C + Mc - M^ c ) = O (Jjj +■■■ = -k,(M% - M^ c )- (4.18) 

If we do not sum up the full corrections, after summing over the polarizations of the outgoing 
massive vector, we would end up with unphysical (non-decoupling) (Ms/mz) 2 and (M x /mz) 2 
corrections [10]. 

4.2 Results for the cross section 

We now turn to evaluate the full 3-body cross section for the process in Eq. (4.1), including 
VIB diagrams. We follow a rather pedagogical approach, starting from the formula for the 
cross section 

* = ^VW* + ~ V%?(*%iW (419) 

being X = [(k\ ■ k 2 ) 2 — M 4 ] 1 / 2 the initial flux. The squared amplitude \M.\ 2 is obtained from 
Eqs. (4.2) and (4.7)-(4.9) by summing over the physical gauge boson polarizations 

E = -^+^- ( 4 - 2 °) 



8=1,2,3 



Z 



Integrating Eq. (4.19) over the three angles that define the position of the plane described 
through the momentum conservation p 1 + p 2 + k = we obtain 

\M\ 2 

vda = ; — -rdx-tdxo , (4.21) 

1024vr 4 ' K J 

where X\ and x 2 parametrize the final energies. In particular, letting s\ = (ki + k 2 ) 2 , we have 

k° = (l-x 2 )ViT/2, (4.22) 
Pi = zivV2, (4-23) 
pi = {l-x x + x 2 )yfc/2, (4.24) 



and we find the following constraints on the phase space 



1 + Xo 1(1 — Xo) 2 rn 2 y 

x_ < xi < x + with x± = — — ± \ (4.25) 

2 y 4 s% 

m| < x 2 <l-2^. (4.26) 



si v s i 
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The integrations of the squared amplitude over the phase space cannot be carried out exactly, 
but two limiting situations are of interest: an expansion in powers of 1/r -C 1, and the case 
with v — > with r generic. The results in the former limit are shown below, while the latter 
case is reported in Appendix B. 

Let us parametrize the cross section as 

I4 fl 



vu 



64tt 2 c^M2 



(Ps + P P ) ■ 



(4.27) 



The partially-inclusive cross section, expanded in the large Ms limit (r 3> 1), is obtained by 
integrating Eq. (4.21) over x\\ neglecting terms vanishing in the mz — > limit we find 



dp s 


4 


dx 2 




dp p 


v 2 


dx 2 


3r 2 



x 2 {l-x 2 f +C(r 



-5\ 



(4.28) 



l + x 2 2 



X-2 



4 11 



+ (1 



X 2 ± 



x 2 



2 - 



r 

2{x 2 



i x + 
log — 



5) 



-3x1 + 14a; 2 



17 x 2 + 34) 



+ £>(r- 5 ). (4.29) 



1 



x 2 



o 

m z 
Ml' 



with x± = 

The fully-inclusive s-wave and p-wave contributions are obtained by further integrating 
over x 2 as prescribed in Eq. (4.26); in the large r limit, we get 

1 



Ps 



pp 



15r 4 

,2 



+ 0(r~ 



v 



180r 2 



+10(15 



7T 



+ 



2 In * - 3 

40(tt 2 - 13) (1059 - llOvr 2 ) 



r / 4 


11> 




60 1 - - + 




\ r 


r 2 y 





+ 



+ 0{r~ 



(4.30) 



(4.31) 



This result shows explicitly the peculiar structure in powers of 1/r that we have estimated 
in Eq. (4.14) using general arguments. As an order-of-magnitude estimate, for velocities 
v ~ 10~ 3 , one expects the 3-body process (in s-wave) to dominate over the 2-body one (in 
p-wave) for values of r < O(10). 

Our results are in disagreement with the ones of Ref. [5], where the helicity suppression 
was removed just by including LO FSR, and where in addition EW corrections proportional 
to M 2 /m?7 were found. Indeed, such terms are present in the physical polarizations sum in 
Eq. (4.20). However, because of gauge invariance in the form of the Ward Identities written 
in Eq. (4.18), they disappear from the final result as they should. On the contrary, we 
find a perfect agreement with the result of Ref. [2] where the photino annihilation process 
77 — > e + e~7 in the v — > limit is analyzed. 



4.3 Effective Field Theory Approach 

Let us now see how the previous results can be interpreted in the effective field theory lan- 
guage. If the mass of the intermediate scalar particle is much larger than the energy scale 
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of the non-relativistic annihilation (E ~ M x ), then it is possible to integrate the scalar out 
and perform an operator expansion in the small parameter 1/r = M^/M|. The effective 
Lagrangian will be given by an infinite series 
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rM 2 



L9l 



+ ..., 



(4.32) 



where (9„ are dimension-n operators. For the theory we considered in Eq. (2.2), the single 
dimension-six operator is 



0, = -\y L \ 2 [xi^x] [hfP L L] , 



(4.33) 



which generates the tree level contributions of Eq. (3.5). The corresponding cross section for 
the XX ff process is 



w(xx -> //)L = 



\Vl\ 



48ttM2 r 2 ' 

A 



(4.34) 



where the usual v 2 -suppression appears. Notice that this result can be recovered from Eq. (3.3) 
in the limit r ^> 1 and corresponds to the estimate in Eq. (3.10). 

As soon as we add the Z emission (for simplicity we consider only the presence of a single 
gauge boson Z), at the LO 0(l/r) in the amplitude only the external leg can radiate a 
gauge boson, making diagrams of the FSR type. We have performed a complete calculation 
using these effective amplitudes. For the total 3-body cross section we find, neglecting terms 
vanishing in the limit mz — > 



va(xX^ffZ)\ 0e = n527r2c ^ 2 



2 V 2 



15-7r 2 + 61n * 

m z 



2 In * 

m z 



. (4.35) 



This expression still bears the v 2 dependence, so that one recovers the result that the EW 
radiation from the external legs cannot remove the p-wave suppression at LO. Notice that 
Eq. (4.35) exhibits the usual single and double logarithmic behavior of infrared origin, and 
that it correctly reproduces the limit for r 3> 1 of Eqs. (4.27) and (4.31). 

From these results it is clear that limiting the analysis to the dimension-six operator in the 
effective theory, which corresponds to work in the LO approximation, misses the right result 
since one could incorrectly conclude that the total cross section itself is p-wave suppressed. As 
shown in the previous section instead, in order to obtain the correct result, namely that the 
cross section receives important s-wave contributions, one needs to consider the diagrams (of 
VIB and FSR type) arising at the next order. Therefore the effect of removing the suppression 
is encoded by operators of dimension higher than six, for example those in C? 8 . 



5 Energy spectra of final stable particles at the interaction point 

The analytical results obtained above have a phenomenological impact for DM indirect 
searches. Indeed, the energy spectra of stable particles produced by DM annihilation, with 
the inclusion of EW bremsstrahlung, can be very different from those commonly obtained by 
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working at the LO. In this section we show our results for the energy spectra at the inter- 
action point, focusing in particular on positrons, antiprotons, photons and neutrinos. The 
propagation of these fluxes of stable particles through the galactic halo will be discussed in 
Section 6. Our analysis is based on the combination of the analytical description of the pri- 
mary annihilation channels with the numerical techniques for subsequent hadronization and 
decay Let us now describe our procedure in more detail. 

As already discussed, we work in the approximation of massless external fermions, under 
which the calculations of the previous sections have been performed. So we consider only the 
case where L = {v^Li^l), for which this is an excellent approximation. In general, channels 
with external fermions of mass rrif would receive other, different contributions proportional to 
(mf/M x ) 2 , in addition to the s-wave contribution from VIB and FSR, as discussed at length 
above. The primary annihilation channels for XX ^> including EW bremsstrahlung, are 

1 = i e L e Li v eL v e L, etell, t\v eL W~ , e~ L v eL W + , e£e£Z, v eL v eL Z} . (5.1) 
The different 3-body channels are simply related by different gauge couplings 

<?{XX -» VeLVeLZ) = \ 2 „ cr(xx -> e£e£Z), (5.2) 

2c 2 

°{XX -> e- L v eL W + ) = a(xx e+u eL W-) = n _ Z2 n 2 ~> 4e^Z)| (5.3) 

<XX^e + L e- Ll ) = (1 4 ^| )2 *{ XX ^etelZ)\ mz ^. (5.4) 

We have written our own Monte Carlo code to generate parton-level events for DM annihila- 
tions into 2- and 3-body final states, in the frame where the total spatial momentum is zero. 
While for the 2 — > 2 processes the final state consists of two back-to-back particles, the 3-body 
final states require particular care because the probability distribution of the momenta of the 
outgoing particles is dictated by the double-differential probability distributions 

1 da(xx ->■ 3-body) 

a(xX ~^ 3— body) dx\&Xi 

where x\ and %2 are related to the energy fractions of the outgoing particles, as in Eqs. (4.22), 
(4.23), (4.24). 

A large number of events (2 x 10 5 ) for each annihilation channel in Eq. (5.1) is generated in 
this way, and then passed through Pythia 8.145 [11] for simulating the subsequent showering, 
hadronization and decay 3 . All unstable particles are requested to decay so that the only final 
particles remaining in the sample are the stable species of the SM. We have performed several 
checks at various levels to assess the reliability of our numerical code. For instance, we have 
found excellent agreement with the results of Ref. [1]. 

3 Pythia 8.1 has been preferred over the predecessor Pythia 6.4 because of the inclusion of the photon 
branchings into fermion-antifermion in the showering process. 
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Figure 3: The spectra dAff/dlnx, as defined in Eq. (5.6), for e + (green), 7 (red), v = {v e -\-v^ + v T ) 
(black) and p (blue), from the annihilation \X ^\^^ v eL^eL with the corresponding weak boson 
emission corrections, for the case M x = 1 TeV, Ms = 4TeV, v = 10~ 3 (solid lines). For comparison, 
we show the spectra (dashed lines) in the LO approximation (see text for details). 



Fitting the numerical results, it is possible to extract the energy distributions of each 
stable particle / 



dM f _ 1 da(xX^f + X) 
dlnx ctq d In x 



f = {e + ,e ,"/,v,v,p,p} , (5.6) 



( f ) ( f) 

where x = E^ etic / M x , E^ etic is the kinetic energy of the particle / (the difference between 
total and kinetic energies is obviously relevant only for the (anti)protons), and the X reminds 
us of the inclusivity in the final state with respect to the particle /. As a normalization, we 
have chosen the tree-level cross section of the 2-body processes 4 

o-q = 0tree(xX -»- e£e^) + a ticc (xx ->■ v eL 9 eL ) . (5.7) 

The plot in Figure 3 shows the resulting dMf/dx for e + ,7, v — (y e + + u T )/3,p for a 
specific, but representative, choice of parameters: M x = 1 TeV, Ms = 4 TeV (corresponding 



4 Another choice for the normalization would be the total cross section a(\X ~^ f + X), which would 
provide the quantity in Eq. (5.6) with a more transparent physical interpretation as the energy spectrum of /. 
However, a(\X —> f + X) is not as easily calculable as the 2-body cross section and it would not be possible 
to compare the energy spectra with and without the s-wave contributions because their total cross sections 
would be different. In any case, the specific choice of the normalization becomes irrelevant when taking ratios 
of spectra, which serve to stress the relevance of the effect we are studying. 
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Primary W 




Figure 4: Ratios between the energy spectrum of the W~ gauge boson in the 3-body annihilation chan- 
ne l XX ~^ ^L v ehW~ for different values of Ms with respect to the same energy spectrum computed 
in the LO approximation (see text for details). 



to r = 16) and v = 1CT 3 ; for comparison, the situation where only the LO term is taken 
into account is also shown (recall that what LO means has been discussed in Sect. 4.1). The 
energy spectra result to be much larger than those obtained in the LO approximation. This 
is a consequence of having a sizeable s-wave annihilation channel opened at the next-to- 
leading order in the 1/r expansion. To better clarify this point we remind that - as already 
discussed in Ref. [1] - the emission of an EW gauge boson opens the hadronic channel, and has 
dramatic consequences on the final state: antiprotons as well as a large number of soft photons, 
positrons and neutrinos from pion decays are produced leading to a huge enhancement in the 
low-energy tail of the energy spectra of final stable particles. The situation described in 
Ref. [1] is obtained here in correspondence of the LO approximation where - as discussed in 
Section 4 - the p-wave term in the 3-body cross section dominates widely over the s-wave one, 
giving corrections factorized with respect to the tree-level annihilation process. Going beyond 
the LO, the opening of a sizeable s-wave contribution results into a harder energy spectrum for 
the primary gauge boson which is entirely converted (after decay and hadronization processes) 
into low energy stable particles, thus leading to a greater enhancement in the low-energy tails 
of their spectra. 

An interesting spectral feature of the gamma rays originates from the inclusion of EW 
corrections (see Figure 3). Indeed, the gamma ray spectrum is the composition of a bump 
in the hard region due to the contribution of hard photons coming from the s-wave in the 
primary annihilation channel XX ~~ e t e Lli anc ^ a huge tail of soft ones originating from 
showering processes and from the hadronization of the W and Z gauge bosons included in 
our analysis. 

The relevance of the effect of removing the suppression is made more manifest by taking the 
ratios between the energy spectra computed at given Ms with respect to those obtained in the 
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Figure 5: Ratios between the energy spectra of final stable particles for different values of M$ with 
respect to those computed in the LO approximation (see text for details). 



LO approximation, as shown in Figure 4 and in Figure 5, for different values of Ms- In Figure 
4 we plot these ratios before hadronization and decays, considering as example the spectrum of 
the primary W~ gauge boson in the 3-body annihilation channel XX ~~ >* ^L v eiW~ ■ The growth 
in the hard region as the s-wave contribution becomes larger is apparent; this prerogative is 
present in all the 3-body channels included in our analysis and listed in Eq. (5.1). In Figure 
5, we show the ratios of the energy spectra of the final stable particles after hadronization 
and decays, with respect to those at LO. The impact of the EW radiation beyond LO can 
result into an enhancement of the energy spectra even by factors O(10 — 100). 

6 Fluxes of final stable particles at detection 

The previous sections have focused on the calculations of the energy spectra of stable SM 
particles at the interaction point, normalized for each DM annihilation event. In this section 
we want to make contact with the phenomenological observables and thus compute the fluxes 
of electrons, positrons, antiprotons, prompt gamma rays and neutrinos that can be measured 
at Earth. We first recall briefly the basics of the computation of such fluxes (see e.g. Ref. [12] 
for a lucid and pedagogical review) and then illustrate the results for the energy spectra found 
in the previous section. 
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6.1 Basics of galactic propagation of stable particles 

Dark Matter distribution in the galactic halo. The DM density profile in the galactic 
halo, p(x), is one of the essential ingredients to determine the normalization of the fluxes of 
cosmic rays that are collected at Earth. N-body numerical simulations performed in the latest 
decades have found different answers for p(r). While recent simulations seem to individuate 
the Einasto profile as the best option, the Navarro-Frenk- White (NFW) profile is stil widely 
used in the literature and the cored Burkert profile (disfavored by simulations) is sometimes 
advocated as a better fit to astronomical observations. These profiles explicitly read 





p{r) 


r s [kpc] 


Ps [GeV/cm 3 ] 


NFW [13] 


Ps r -(l+ r - 

r V r sj 


-2 

1 


24.42 


0.184 


Einasto [14] 


^ exp [-637 [(£) 


0.17 ~T 
- 1 


28.44 


0.033 


Burkert [15] 


Ps 

(l + r/r s )(l + (r/r s ) 2 ) 


12.7 


0.712 



where, in order to fix the parameters r s and p s at their precise values, one imposes the 
constraints that p & (the value of the DM density at the location of the solar system) = 
0.3 GeV/cm 3 and that the total DM mass contained in the Milky Way reproduces observations 
(see Ref. [16]). They differ most at the Galactic Center (GC): NFW is peaked as r~ l while 
Burkert is constant in the inner 1 kpc. They are instead similar around the location of the 
solar system, due also to the p & constraint. As long as a convergent determination of the 
actual DM profile is not reached, it is sensible to have at disposal the whole range of these 
possible choices when computing Dark Matter signals in the Milky Way. In other words, 
the ignorance on the actual DM profile constitutes a (currently) irreducible astrophysical 
uncertainty for the predicted fluxes. 

Charged particles (electrons, positrons, antiprotons). The e~, e + and p produced in 
any given point of the halo propagate immersed in the turbulent galactic magnetic field. The 
field consists of random inhomogeneities that act as scattering centers for charged particles, so 
that their journey can effectively be described as a diffusion process from an extended source 
(the DM halo) to some final given point (the location of the Earth, in the case of interest). 
The number density rif(x,E) per unit energy E of the cosmic ray species / (= e + ,e~,p) in 
any given point x evolves according to a diffusion-loss equation [12] 

The first term accounts for diffusion, with a coefficient conventionally parameterized as 
)C(E) = Kq{EJ GeV) 5 . The second term describes energy losses: the coefficient b is position- 
dependent since the intensity of the magnetic field (which determines losses due to synchrotron 
radiation) and the distribution of the photon field (which determines losses due to inverse 
Compton scattering) vary across the galactic halo. The third term deals with convection 
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while the last term accounts for nuclear spallations, that occur with rate V in the disk of 
thickness h ~ 100 pc. The source, DM annihilations, is denoted by Q. The different pro- 
cesses described above have a different importance depending on the particle species: the 
journey of electrons and positrons is primarily affected by synchrotron radiation and inverse 
Compton energy losses, while for antiprotons these losses are negligible and convection and 
spallation dominate. 

Eq. (6.1) is usually solved numerically in a diffusive region with the shape of a solid 
flat cylinder that sandwiches the galactic plane, with height 2L in the z direction and radius 
R = 20 kpc in the r direction. The location of the solar system corresponds to x Q = (r , z@) = 
(8.33 kpc, 0). Boundary conditions are imposed such that the number density rif vanishes on 
the surface of the cylinder, outside of which the charged cosmic rays freely propagate and 
escape. The values of the propagation parameters 5, K , V conv and L are deduced from 
a variety of (ordinary) cosmic ray data and modelizations. It is customary to adopt the 
following sets, denoted with MIN, MED and MAX because they are found to minimize or 
maximize the final fluxes 





Electrons or positrons 




Antiprotons 
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5 /C [kpc 2 /Myr] 


5 


/C [kpc 2 /Myr] V conv [km/s] 


L [kpc] 


MIN 


0.55 0.00595 


0.85 
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1 


MED 
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0.0112 12 


4 


MAX 


0.46 0.0765 


0.46 
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As long as independent measurements do not allow to pin down more precisely the values of 
these parameters, the scatter among such different sets constitute an additional astrophysical 
uncertainty on the predicted DM fluxes, this time due to the propagation process. 

The solution of Eq. (6.1) allows to compute the phenomenological quantity in which we 
are interested: the flux of cosmic rays received at Earth d&f/dE = VfUf/An (where Vf is 
the velocity of species /, equal to c for e but possibly different for mildly-relativistic p). It 
turns out that, both for e and for p, the flux can be conveniently expressed as a convolution 
of the spectra at the interaction point with some universal functions that encapsulate the 
astrophysics of the 'production and propagation' process. More precisely, for e ± one has 



(E,x e 



v P ± 



dE v '~ 0j Airb(E,x G )2 \M X/ 

where dN e ±/dE are the spectra at the annihilation point and I(E, E s ,x & ) are (generalized) 
halo functions which are independent of the particle physics model: there is such a function 
for each choice of DM distribution profile and choice of e ± propagation parameters. We are 
following here the formalism discussed in Ref. [16], which allows in particular to take into 
account the spatial dependence of the energy loss coefficient b for discussed above. We 
refer to Ref. [16] for all details, including an explicit form of b and of the I functions, and 
further references. Similarly, for p one has 

2 



Pe 



(av) 



M 



dE 



(E s ) I(E, E s , x e 



(6.2) 
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where it is now the function R(E) which contains the astrophysics: again, there is such a 
function for each choice of DM distribution profile and the choice of p propagation parameters. 

Neutral particles (photons, neutrinos). Neutral messengers produced by DM annihi- 
lation in any given point of the DM halo travel along a straight line to the Earth. Since 
absorption in the Galaxy is negligible, the flux from a given direction is the result of the 
contribution from all the Dark Matter intervening along the line of sight. The integrated flux 
of gamma rays or neutrinos over a region AQ, corresponding e.g. to the window of observation 
or the resolution of the telescope, is given by 

d 4¥(E) = ?-(^) 2 J^(av) d -^, withJ^- 1 -/ [ ^ (rM)) V,(6.4) 
dE 1 ; An2\M x J dE AO J An J ho s r \ Po ) A 1 

where dN ltV /dE denotes as usual the spectra at the annihilation point and the average J 
factor contains the integral along the line of sight (l.o.s.). Here the coordinate r, centered on 
the GC, reads r(s, 6) = (r + s 2 — 2 r s cos #) 1//2 : s runs along the l.o.s. and is the aperture 
angle between the direction of the l.o.s. and the axis connecting the Earth to the GC. For a 
fixed window AQ, the value of J can span orders of magnitude depending on the choice of the 
DM profile, especially if the window is small and close to the region where the profiles differ 
most, i.e. the GC. We refer to Ref. [16] for some explicit values of J in selected windows. 



6.2 Results 

In Figure 6 we show the fluxes of charged cosmic rays for the representative choice of model 
parameters already used in Section 5 (namely, M x = 1 TeV and Ms = 4 TeV). We do not 
adopt here a specific value for the annihilation cross section (<jv), which obviously enters as 
a normalizing constant in Eqs. (6.2) and (6.3): in the context of the toy model that we are 
considering, its value is very small, if one requires the cosmological relic abundance of this 
DM candidate to be fixed by the thermal freeze-out mechanism. This is not surprising and it 
is actually the typical case for Bino-like Dark Matter in supersymmetry. Notice that if one 
adopts the normalization with the 2-body annihilations, as done in Sect. 5, the energy spectra 
and the cross section appearing in the expressions for the fluxes should be replaced by the 
analogous quantities in Eqs. (5.6) and (5.7). What we are more interested in is verifying that 
the enhancement of the fluxes remains significant on the phenomenological observables (the 
spectra after propagation) and distinct from the normalization issues due to the propagation 
itself. This is indeed the Figure 6 shows. In these plots, the solid lines represent 

the fluxes computed with 'NFW as the choice of DM profile and with 'MED' propagation 
parameters. The shaded bands show the variance of this prediction that one obtains by 
making other choices. The bands span quite a large area since both the DM profile and the 
propagation parameters are variated simultaneously. 

The fluxes of e and p are about one order of magnitude higher than those computed in 
the LO approximation, consistently with what expected from input fluxes, and emerge quite 
clearly from the uncertainty bands, especially for electrons or positrons at high energies. Small 
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Electrons or Positrons 



Antiprotons 




differences in the shapes of the spectra, that were marginally visible in the input spectra, see 
Figure 3, are essentially washed out for electrons and positrons (since the energy losses tend 
to smooth any spectral feature) but remain somewhat discernible in the antiproton spectrum 
(since propagation does not significantly reshuffle energies for this species). 

In the case of neutral particles (gamma rays and neutrinos), the fluxes 'at detection' are 
easily computed with the use of Eq. (6.4). They simply correspond to a re-normalization of 
the input fluxes, for a given choice of the observational window and the DM profile (which 
fixes the J factor), so that we do not plot them explicitly. Any peculiar spectral feature 
possibly introduced by contributions beyond the LO approximation, e.g. the high-energy 
bump discussed above in the gamma ray spectrum, would of course be conserved. 



7 Summary and conclusions 

We have investigated the relevance of the EW corrections in theories where the cross section 
for DM annihilation into 2-body final states is suppressed. A Majorana DM annihilating into 
two light SM fermions is one such case. We have worked for simplicity with a model where the 
DM is a Majorana fermion of mass M x and a SM singlet, which annihilates into SM fermions 
through the exchange of a heavy scalar doublet of mass Ms, and carried out an expansion in 
l/r = {M x /M s f « 1. 

Let us summarize our main results: 
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• at the lowest order (1/r in the amplitude) the radiation of EW gauge bosons is not able 
to remove the helicity suppression and the process stays in the p-wave (see Eq. (3.10) 
for an estimate and Eq. (3.3) for the precise result); 

• an efficient removal of the suppression, opening up a potentially large s-wave, is achieved 
by including EW radiation at the next-to- leading order (1/r 2 in the amplitude), which 
comes from both FSR and VIB diagrams (see Eq. (4.14) for an estimate and Eqs. (4.30)- 
(4.31) for the precise result); 

• the resulting energy spectra of stable particles, in the annihilation region, get substan- 
tially enhanced by this effect by factors O(10 — 100) (see Figure 5). 

Furthermore, such an effect does not get spoiled by galactic propagation and crucially affects 
the predictions for fluxes to be measured at Earth (see Figure 6). 

We have also interpreted our findings in the language of effective field theory and pointed 
out that the effect of opening up the s-wave is missed by dimension-six operators and only 
catched by higher-dimensional operators. This is an example where the naive dimensional 
power counting fails to assess the relative importance of the operators in the expansion, as 
far as the EW radiation is concerned. 

Our results have a wider generality than the specific model we have considered. Reliable 
computations of energy spectra of stable particles and predictions for their fluxes at Earth - 
the key observable for DM indirect searches - cannot prescind from including the effects of 
EW radiation. 
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A The Dirac case 

We list here the analogous results of Sections 3 and 4 for the case of Dirac Dark Matter. 
The 2-body annihilation cross section into a pair of massless left-handed fermions is 



va = a + bv 2 + 0(v 4 ) , 



(A.l) 



where 



a = 



32tt(1 + r) 2 M 2 ' 




b = 



y L \ 4 (r 2 -3r- 1) 
96ttM 2 (1 + r) 4 



(A.2) 
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For the process in Eq. (4.1), the exchanged diagrams of Figure 2 are absent and the amplitude 
can be written as 

ig\y L \ 2 (l-2s 2 w ) 



iM-e* = 



[M A + M b + M c ] , 



(A.3) 



and by using the Fierz transformation (4.4) the three terms analog of Eqs. (4.7)-(4.9) become 
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Then the calculation of the cross sections proceeds as described in Section 4.2 and we choose 
the parametrization in terms of p s and p p as in Eq. (4.27). In particular for the partially- 
inclusive cross section in the large Ms limit we find, neglecting terms vanishing in the mz — > 
limit 
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For completeness, we also evaluate the expressions for the fully-inclusive s-wave and p-wave 
contributions 
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B 3-body cross section in the v — > limit 



We report here the results for the cross section of the 3-body process XX ~~ ffZ, in the limit 
v — > 0, therefore retaining only the part of the process proceeding through the s-wave. We 
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do not expand in powers of 1/r, so the following results are valid for any value of r > 1. The 
cross section is parametrized as in Eq. (4.27) 
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The partially-inclusive contribution to the differential cross section is 
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with y = y (1 — x 2 ) 2 — ^gf . The fully-inclusive cross section is obtained by integrating the 
previous expressions over the kinematical domain in Eq. (4.26). Neglecting terms vanishing 
in the limit mz — > we find 
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being Li 2 (z) = YlkLi ^V^ 2 the usual dilogarithm. 
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